home *** CD-ROM | disk | FTP | other *** search
/ Space & Astronomy / Space and Astronomy (October 1993).iso / mac / programs / space / DE118I.ZIP / ADAMS4.C < prev    next >
C/C++ Source or Header  |  1993-02-06  |  17KB  |  717 lines
  1. /* Adams-Bashforth-Moulton integration formulas.
  2.  *
  3.  * Reference:
  4.  * Shampine, L. F. and M. K. Gordon, _Computer Solution of
  5.  * Ordinary Differential Equations_, W. H. Freeman, 1975.
  6.  *
  7.  * Program by Steve Moshier.
  8.  */
  9.  
  10. #ifndef NOINCS
  11. #include "mconf.h"
  12. #include "prec.h"
  13. #include "int.h"
  14. #include "const.h"
  15. #endif
  16.  
  17. /* Divided differences */
  18. #define N 19
  19.  
  20.  
  21. /*Predictor coefficients:
  22.  1.0
  23.  1.0 / 2.0
  24.  5.0 / 12.0
  25.  3.0 / 8.0
  26.  251.0 / 720.0
  27.  95.0 / 288.0
  28.  19087.0 / 60480.0
  29.  5257.0 / 17280.0
  30.  1070017.0 / 3628800.0
  31.  25713.0 / 89600.0
  32.  26842253.0 / 95800320.0
  33.  4777223.0 / 17418240.0
  34.  703604254357.0 / 2615348736000.0
  35.  106364763817.0 / 402361344000.0
  36.  1166309819657.0 / 4483454976000.0
  37.  2.5221445e7 / 9.8402304e7
  38.  8.092989203533249e15 / 3.201186852864e16
  39.  8.5455477715379e13 / 3.4237292544e14
  40.  1.2600467236042756559e19 / 5.109094217170944e19
  41. */
  42.  
  43. /* Corrector coefficients:
  44.    1.0
  45.   -1.0 / 2.0
  46.   -1.0 / 12.0
  47.   -1.0 / 24.0
  48.   -19.0 / 720.0
  49.   -3.0 / 160.0
  50.   -863.0 / 60480.0
  51.   -275.0 / 24192.0
  52.   -33953.0 / 3628800.0
  53.   -8183.0 / 1036800.0
  54.   -3250433.0 / 479001600.0
  55.   -4671.0 / 788480.0
  56.   -13695779093.0 / 2615348736000.0
  57.   -2224234463.0 / 475517952000.0
  58.   -132282840127.0 / 31384184832000.0
  59.   -2639651053.0 / 689762304000.0
  60.   1.11956703448001e14 / 3.201186852864e16
  61.   5.0188465e7 / 1.5613165568e10
  62.   2.334028946344463e15 / 7.86014494949376e17
  63. */
  64. #if LDOUBLE
  65. #if IBMPC
  66. static short precof[] = {
  67. 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD
  68. 0x0000,0x0000,0x0000,0x8000,0x3ffe, XPD
  69. 0x5555,0x5555,0x5555,0xd555,0x3ffd, XPD
  70. 0x0000,0x0000,0x0000,0xc000,0x3ffd, XPD
  71. 0xd27d,0x7d27,0x27d2,0xb27d,0x3ffd, XPD
  72. 0x38e4,0xe38e,0x8e38,0xa8e3,0x3ffd, XPD
  73. 0x5440,0xea99,0x43fe,0xa195,0x3ffd, XPD
  74. 0x3519,0x6dfc,0x518a,0x9bc3,0x3ffd, XPD
  75. 0xb126,0x0fb1,0xf045,0x96f8,0x3ffd, XPD
  76. 0x80bb,0x54d8,0x721a,0x92ee,0x3ffd, XPD
  77. 0x6246,0xcff2,0x02c2,0x8f75,0x3ffd, XPD
  78. 0xac15,0xb603,0x8869,0x8c6c,0x3ffd, XPD
  79. 0xdc52,0x2596,0x2625,0x89be,0x3ffd, XPD
  80. 0x2820,0xc6b2,0x0f57,0x8759,0x3ffd, XPD
  81. 0xf898,0xbc14,0x9903,0x8530,0x3ffd, XPD
  82. 0x4d10,0xe1e8,0xff92,0x833a,0x3ffd, XPD
  83. 0x883c,0x772e,0x97fc,0x8170,0x3ffd, XPD
  84. 0xe995,0xbf9f,0x86c4,0xff96,0x3ffc, XPD
  85. 0x78b9,0x72bf,0x1a81,0xfc8c,0x3ffc, XPD
  86. };
  87. static short corcof[] = {
  88. 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD
  89. 0x0000,0x0000,0x0000,0x8000,0xbffe, XPD
  90. 0xaaab,0xaaaa,0xaaaa,0xaaaa,0xbffb, XPD
  91. 0xaaab,0xaaaa,0xaaaa,0xaaaa,0xbffa, XPD
  92. 0xd82e,0x2d82,0x82d8,0xd82d,0xbff9, XPD
  93. 0x999a,0x9999,0x9999,0x9999,0xbff9, XPD
  94. 0x9474,0x1e9c,0x473f,0xe9c9,0xbff8, XPD
  95. 0xe4e9,0x93a3,0x4e8f,0xba3e,0xbff8, XPD
  96. 0x7e46,0xc950,0x28ab,0x994c,0xbff8, XPD
  97. 0x0d67,0x5b26,0xc557,0x814f,0xbff8, XPD
  98. 0x9d4e,0x3987,0xd5e1,0xde5b,0xbff7, XPD
  99. 0x8c4d,0x7bad,0x9646,0xc21e,0xbff7, XPD
  100. 0xf0ac,0x1b33,0x9124,0xab98,0xbff7, XPD
  101. 0x0c7c,0xb92d,0xb357,0x9945,0xbff7, XPD
  102. 0xe201,0xa74b,0x9502,0x8a1d,0xbff7, XPD
  103. 0xc3dc,0x1655,0xb86d,0xfacc,0xbff6, XPD
  104. 0x6a55,0x5ce2,0xcb35,0xe533,0x3ff6, XPD
  105. 0xb891,0xaf49,0x4d0b,0xd2aa,0x3ff6, XPD
  106. 0x370d,0x381c,0x10d3,0xc29b,0x3ff6, XPD
  107. };
  108. #endif
  109. #if MIEEE
  110. static long precof[3*19] = {
  111. 0x3fff0000,0x80000000,0x00000000,
  112. 0x3ffe0000,0x80000000,0x00000000,
  113. 0x3ffd0000,0xd5555555,0x55555555,
  114. 0x3ffd0000,0xc0000000,0x00000000,
  115. 0x3ffd0000,0xb27d27d2,0x7d27d27d,
  116. 0x3ffd0000,0xa8e38e38,0xe38e38e4,
  117. 0x3ffd0000,0xa19543fe,0xea995440,
  118. 0x3ffd0000,0x9bc3518a,0x6dfc3519,
  119. 0x3ffd0000,0x96f8f045,0x0fb1b126,
  120. 0x3ffd0000,0x92ee721a,0x54d880bb,
  121. 0x3ffd0000,0x8f7502c2,0xcff26246,
  122. 0x3ffd0000,0x8c6c8869,0xb603ac15,
  123. 0x3ffd0000,0x89be2625,0x2596dc52,
  124. 0x3ffd0000,0x87590f57,0xc6b22820,
  125. 0x3ffd0000,0x85309903,0xbc14f898,
  126. 0x3ffd0000,0x833aff92,0xe1e84d10,
  127. 0x3ffd0000,0x817097fc,0x772e883c,
  128. 0x3ffc0000,0xff9686c4,0xbf9fe995,
  129. 0x3ffc0000,0xfc8c1a81,0x72bf78b9,
  130. };
  131. static long corcof[3*19] = {
  132. 0x3fff0000,0x80000000,0x00000000,
  133. 0xbffe0000,0x80000000,0x00000000,
  134. 0xbffb0000,0xaaaaaaaa,0xaaaaaaab,
  135. 0xbffa0000,0xaaaaaaaa,0xaaaaaaab,
  136. 0xbff90000,0xd82d82d8,0x2d82d82e,
  137. 0xbff90000,0x99999999,0x9999999a,
  138. 0xbff80000,0xe9c9473f,0x1e9c9474,
  139. 0xbff80000,0xba3e4e8f,0x93a3e4e9,
  140. 0xbff80000,0x994c28ab,0xc9507e46,
  141. 0xbff80000,0x814fc557,0x5b260d67,
  142. 0xbff70000,0xde5bd5e1,0x39879d4e,
  143. 0xbff70000,0xc21e9646,0x7bad8c4d,
  144. 0xbff70000,0xab989124,0x1b33f0ac,
  145. 0xbff70000,0x9945b357,0xb92d0c7c,
  146. 0xbff70000,0x8a1d9502,0xa74be201,
  147. 0xbff60000,0xfaccb86d,0x1655c3dc,
  148. 0x3ff60000,0xe533cb35,0x5ce26a55,
  149. 0x3ff60000,0xd2aa4d0b,0xaf49b891,
  150. 0x3ff60000,0xc29b10d3,0x381c370d,
  151. };
  152. #endif /* ldouble MIEEE */
  153. #else /* LDOUBLE */
  154. #if UNK
  155. static DOUBLE precof[19] = {
  156.  1.00000000000000000000E0,
  157.  5.00000000000000000000E-1,
  158.  4.16666666666666666667E-1,
  159.  3.75000000000000000000E-1,
  160.  3.48611111111111111111E-1,
  161.  3.29861111111111111111E-1,
  162.  3.15591931216931216931E-1,
  163.  3.04224537037037037037E-1,
  164.  2.94868000440917107584E-1,
  165.  2.86975446428571428571E-1,
  166.  2.80189596443936721714E-1,
  167.  2.74265540031599059377E-1,
  168.  2.69028846773648774310E-1,
  169.  2.64351348366606509794E-1,
  170.  2.60136396127601036937E-1,
  171.  2.56309496574389152514E-1,
  172.  2.52812146729039234860E-1,
  173.  2.49597650297715667409E-1,
  174.  2.46628202582257457797E-1,
  175. };
  176. static DOUBLE corcof[19] = {
  177.  1.00000000000000000000E0,
  178. -5.00000000000000000000E-1,
  179. -8.33333333333333333333E-2,
  180. -4.16666666666666666667E-2,
  181. -2.63888888888888888889E-2,
  182. -1.87500000000000000000E-2,
  183. -1.42691798941798941799E-2,
  184. -1.13673941798941798942E-2,
  185. -9.35653659611992945326E-3,
  186. -7.89255401234567901235E-3,
  187. -6.78584998463470685693E-3,
  188. -5.92405641233766233766E-3,
  189. -5.23669325795028506669E-3,
  190. -4.67749840704226451581E-3,
  191. -4.21495223900547285688E-3,
  192. -3.82689955321188442330E-3,
  193.  3.49734984534991765411E-3,
  194.  3.21449643132356745146E-3,
  195.  2.96944771545820961119E-3,
  196. };
  197. #endif /* double UNK */
  198. #if DEC
  199. static short precof[19*4] = {
  200. 0040200,0000000,0000000,0000000,
  201. 0040000,0000000,0000000,0000000,
  202. 0037725,0052525,0052525,0052525,
  203. 0037700,0000000,0000000,0000000,
  204. 0037662,0076447,0151175,0023722,
  205. 0037650,0161616,0034343,0107071,
  206. 0037641,0112503,0177352,0114524,
  207. 0037633,0141521,0105155,0176065,
  208. 0037626,0174360,0042417,0130661,
  209. 0037622,0167162,0015124,0154201,
  210. 0037617,0072402,0141317,0171142,
  211. 0037614,0066210,0064666,0001654,
  212. 0037611,0137046,0022445,0113334,
  213. 0037607,0054417,0053706,0131050,
  214. 0037605,0030231,0001674,0012371,
  215. 0037603,0035377,0111341,0164115,
  216. 0037601,0070227,0176167,0027210,
  217. 0037577,0113206,0142277,0117752,
  218. 0037574,0106032,0100562,0137571,
  219. };
  220. static short corcof[19*4] = {
  221. 0040200,0000000,0000000,0000000,
  222. 0140000,0000000,0000000,0000000,
  223. 0137252,0125252,0125252,0125253,
  224. 0137052,0125252,0125252,0125253,
  225. 0136730,0026602,0154055,0101330,
  226. 0136631,0114631,0114631,0114632,
  227. 0136551,0144507,0037436,0116224,
  228. 0136472,0037116,0107623,0121745,
  229. 0136431,0046050,0125711,0050176,
  230. 0136401,0047705,0053533,0023015,
  231. 0136336,0055725,0160471,0103635,
  232. 0136302,0017226,0043173,0126614,
  233. 0136253,0114221,0022033,0031761,
  234. 0136231,0042663,0053671,0026414,
  235. 0136212,0016625,0001247,0045742,
  236. 0136172,0146270,0066426,0052704,
  237. 0036145,0031713,0032534,0161152,
  238. 0036122,0125115,0005657,0044671,
  239. 0036102,0115420,0151470,0016067,
  240. };
  241. #endif /* double DEC */
  242. #if IBMPC
  243. static short precof[19*4] = {
  244. 0x0000,0x0000,0x0000,0x3ff0,
  245. 0x0000,0x0000,0x0000,0x3fe0,
  246. 0xaaab,0xaaaa,0xaaaa,0x3fda,
  247. 0x0000,0x0000,0x0000,0x3fd8,
  248. 0xa4fa,0xfa4f,0x4fa4,0x3fd6,
  249. 0x71c7,0xc71c,0x1c71,0x3fd5,
  250. 0x532b,0x7fdd,0x32a8,0x3fd4,
  251. 0xbf87,0x314d,0x786a,0x3fd3,
  252. 0xf636,0x08a1,0xdf1e,0x3fd2,
  253. 0x9b10,0x434a,0x5dce,0x3fd2,
  254. 0xfe4c,0x5859,0xeea0,0x3fd1,
  255. 0xc076,0x0d36,0x8d91,0x3fd1,
  256. 0xb2dc,0xc4a4,0x37c4,0x3fd1,
  257. 0xd645,0xeaf8,0xeb21,0x3fd0,
  258. 0x829f,0x2077,0xa613,0x3fd0,
  259. 0x3d0a,0xf25c,0x675f,0x3fd0,
  260. 0xe5d1,0xff8e,0x2e12,0x3fd0,
  261. 0xf3fd,0xd897,0xf2d0,0x3fcf,
  262. 0x57ef,0x502e,0x9183,0x3fcf,
  263. };
  264. static short corcof[19*4] = {
  265. 0x0000,0x0000,0x0000,0x3ff0,
  266. 0x0000,0x0000,0x0000,0xbfe0,
  267. 0x5555,0x5555,0x5555,0xbfb5,
  268. 0x5555,0x5555,0x5555,0xbfa5,
  269. 0xb05b,0x5b05,0x05b0,0xbf9b,
  270. 0x3333,0x3333,0x3333,0xbf93,
  271. 0xd393,0xe7e3,0x3928,0xbf8d,
  272. 0x747d,0xd1f2,0x47c9,0xbf87,
  273. 0x2a10,0x1579,0x2985,0xbf83,
  274. 0x64c2,0xaaeb,0x29f8,0xbf80,
  275. 0x30f4,0xbc27,0xcb7a,0xbf7b,
  276. 0x75b2,0xc8cf,0x43d2,0xbf78,
  277. 0x667e,0x2483,0x7312,0xbf75,
  278. 0x25a2,0x6af7,0x28b6,0xbf73,
  279. 0xe97c,0xa054,0x43b2,0xbf71,
  280. 0xcab8,0x0da2,0x5997,0xbf6f,
  281. 0x9c4d,0x66ab,0xa679,0x3f6c,
  282. 0xe937,0xa175,0x5549,0x3f6a,
  283. 0x0387,0x1a67,0x5362,0x3f68,
  284. };
  285. #endif /* double IBMPC */
  286. #if MIEEE
  287. static short precof[19*4] = {
  288. 0x3ff0,0x0000,0x0000,0x0000,
  289. 0x3fe0,0x0000,0x0000,0x0000,
  290. 0x3fda,0xaaaa,0xaaaa,0xaaab,
  291. 0x3fd8,0x0000,0x0000,0x0000,
  292. 0x3fd6,0x4fa4,0xfa4f,0xa4fa,
  293. 0x3fd5,0x1c71,0xc71c,0x71c7,
  294. 0x3fd4,0x32a8,0x7fdd,0x532b,
  295. 0x3fd3,0x786a,0x314d,0xbf87,
  296. 0x3fd2,0xdf1e,0x08a1,0xf636,
  297. 0x3fd2,0x5dce,0x434a,0x9b10,
  298. 0x3fd1,0xeea0,0x5859,0xfe4c,
  299. 0x3fd1,0x8d91,0x0d36,0xc076,
  300. 0x3fd1,0x37c4,0xc4a4,0xb2dc,
  301. 0x3fd0,0xeb21,0xeaf8,0xd645,
  302. 0x3fd0,0xa613,0x2077,0x829f,
  303. 0x3fd0,0x675f,0xf25c,0x3d0a,
  304. 0x3fd0,0x2e12,0xff8e,0xe5d1,
  305. 0x3fcf,0xf2d0,0xd897,0xf3fd,
  306. 0x3fcf,0x9183,0x502e,0x57ef,
  307. };
  308. static short corcof[19*4] = {
  309. 0x3ff0,0x0000,0x0000,0x0000,
  310. 0xbfe0,0x0000,0x0000,0x0000,
  311. 0xbfb5,0x5555,0x5555,0x5555,
  312. 0xbfa5,0x5555,0x5555,0x5555,
  313. 0xbf9b,0x05b0,0x5b05,0xb05b,
  314. 0xbf93,0x3333,0x3333,0x3333,
  315. 0xbf8d,0x3928,0xe7e3,0xd393,
  316. 0xbf87,0x47c9,0xd1f2,0x747d,
  317. 0xbf83,0x2985,0x1579,0x2a10,
  318. 0xbf80,0x29f8,0xaaeb,0x64c2,
  319. 0xbf7b,0xcb7a,0xbc27,0x30f4,
  320. 0xbf78,0x43d2,0xc8cf,0x75b2,
  321. 0xbf75,0x7312,0x2483,0x667e,
  322. 0xbf73,0x28b6,0x6af7,0x25a2,
  323. 0xbf71,0x43b2,0xa054,0xe97c,
  324. 0xbf6f,0x5997,0x0da2,0xcab8,
  325. 0x3f6c,0xa679,0x66ab,0x9c4d,
  326. 0x3f6a,0x5549,0xa175,0xe937,
  327. 0x3f68,0x5362,0x1a67,0x0387,
  328. };
  329. #endif /* double MIEEE */
  330. #endif /* not LDOUBLE */
  331.  
  332. /* Compute zeroth through kth backward differences
  333.  * of the data in the input array
  334.  */
  335. divdif(vec , k, diffn)
  336. DOUBLE vec[]; /* input array of k+1 data items */
  337. DOUBLE *diffn; /* output array of ith differences */
  338. int k;
  339. {
  340. DOUBLE diftbl[N];
  341. DOUBLE *p, *q;
  342. DOUBLE y;
  343. int i, o;
  344.  
  345. /* Copy the given data (zeroth difference) into temp array
  346.  */
  347. p = diftbl;
  348. q = vec;
  349. for( i=0; i<=k; i++ )
  350.     *p++ = *q++;
  351.  
  352. /* On the first outer loop, k-1 first differences are calculated.
  353.  * These overwrite the original data in the temp array.
  354.  */
  355. o = k;
  356. for( o=k; o>0; o-- )
  357.     {
  358.     p = diftbl;
  359.     q = p;
  360.     for( i=0; i<o; i++ )
  361.         {
  362.         y = *p++;
  363.         *q++ = *p - y;
  364.         }
  365.     *diffn++ = *p; /* copy out the last (undifferenced) item */
  366. #if DEBUG
  367.     printf( "%.5e ", *p );
  368. #endif
  369.     }
  370. #if DEBUG
  371.     printf( "%.5e\n", *(q-1) );
  372. #endif
  373. *diffn++ = *(q-1);
  374. }
  375.  
  376.  
  377. /* Update array of differences, given new data value.
  378.  * diffn is an array of k+1 differences, starting with the
  379.  * zeroth difference (the previous original data value).
  380.  */
  381. dupdate( diffn, k, f )
  382. register DOUBLE *diffn;  /* input and output array of differences */
  383. int k; /* max order of differences */
  384. DOUBLE f; /* new data point (zeroth difference) */
  385. {
  386. DOUBLE new, old;
  387. int i;
  388.  
  389. new = f;
  390. for( i=0; i<k; i++ )
  391.     {
  392.     old = *diffn;
  393.     *diffn++ = new;
  394. #if DEBUG
  395.     printf( "%.5e ", new );
  396. #endif
  397.     new = new - old;
  398.     }
  399. #if DEBUG
  400.     printf( "%.5e\n", new );
  401. #endif
  402. *diffn++ = new;
  403. }
  404.  
  405.  
  406.  
  407.  
  408. /* Evaluate the interpolating polynomial
  409.  *
  410.  *              (x - x )
  411.  *                    n    1
  412.  * P(x) = f  +  --------  D f   +  ...
  413.  *         n       h         n
  414.  *
  415.  *     (x - x )(x - x   )...(x - x     )
  416.  *           n       n-1          n+2-k    k-1
  417.  *  +  ---------------------------------  D    f
  418.  *                   k-1                        n
  419.  *                  h     (k-1)!
  420.  *
  421.  *
  422.  *         j
  423.  *  where D denotes the jth backward difference, see dupdate(), and
  424.  *
  425.  *  f   =   f( x , y(x ) )  is the interpolated derivative y'(x ) .
  426.  *   n          n     n                                        n
  427.  *
  428.  * The subroutine argument t is linearly scaled so that t = 1.0
  429.  * will evaluate the polynomial at x = x_n + h,
  430.  * t = 0.0 corresponds to x = x_n, etc.
  431.  */
  432. DOUBLE difpol( diffn, k, t )
  433. DOUBLE *diffn;
  434. int k; /* differences go up to order k-1 */
  435. DOUBLE t; /* scaled argument */
  436. {
  437. DOUBLE f, fac, s, u;
  438. int i;
  439.  
  440. f = *diffn++; /* the zeroth difference = nth data point */
  441. u = One;
  442. /*s = x/h - n;*/
  443. s = One;    /* to evaluate the polynomial at x = xn + h */
  444. fac = One;
  445. for( i=1; i<k; i++ )
  446.     {
  447.     if( s == 0 )
  448.         break;
  449.     u *= s / fac;
  450.     f += u  *  (*diffn++);
  451.     fac += One;
  452.     s += One;
  453.     }
  454. return( f );
  455. }
  456.  
  457.  
  458. /* Integrate the interpolating polynomial from x[n] to x[n+1]
  459.  * to obtain the change in the integrated function from y[n] to y[n+1]
  460.  * given by:
  461.  *
  462.  *                     k-1
  463.  *                      -       i
  464.  *  y    =  y   +   h   >   c  D  f    .
  465.  *   n+1     n          -    i     n
  466.  *                     i=0
  467.  *
  468.  *  This subroutine returns the summation term, not multiplied by h.
  469.  *  The coefficients c_i of the integration formula are either
  470.  *  precof[] or corcof[], given above.
  471.  */
  472. DOUBLE intpol( diffn, coeffs, k )
  473. DOUBLE *diffn; /* array of backward differences */
  474. DOUBLE *coeffs; /* coefficients of integration formula */
  475. int k; /* differences used go up to order k-1 */
  476. {
  477. DOUBLE s;
  478. int i;
  479.  
  480. s = Zero;
  481. coeffs += k;
  482. diffn += k;
  483. for( i=0; i<k; i++ )
  484.     {
  485.     s += (*--coeffs) * (*--diffn);
  486.     }
  487. return( s );
  488. }
  489.  
  490.  
  491.  
  492. /* Copy array of n elements from p to q.
  493.  */
  494. vcopy( p, q, n )
  495. register DOUBLE *p, *q;
  496. register int n;
  497. {
  498. do
  499.     *q++ = *p++;
  500. while( --n );
  501. }
  502.  
  503.  
  504. /* Adams initialization program.
  505.  */
  506.  
  507. /* Addresses within the work array */
  508. DOUBLE *dv;
  509. DOUBLE *dvp;
  510. DOUBLE *vp;
  511. DOUBLE *yp;
  512. DOUBLE *vn;
  513. DOUBLE *delta;
  514. DOUBLE *sdelta;
  515. DOUBLE *adamy0;
  516.  
  517. DOUBLE ccor;
  518. DOUBLE hstep; /* step size (constant) */
  519. int order; /* Order of the prediction formula */
  520. int ordp1;
  521. int asiz;
  522. int dsiz;
  523. int ineq;
  524. long jstep; /* counts steps taken */
  525.  
  526.  
  527. /* Initialize pointers in work array, compute derivatives at
  528.  * initial position, and start the difference tables.
  529.  * neq >= nequat in adstep() below.  If neq > nequat, unused space
  530.  * is left for extra variables that are not actually integrated.
  531.  */
  532. adstart( h, yn, work, neq, ord, t )
  533. DOUBLE h;
  534. DOUBLE yn[], work[];
  535. int neq, ord;
  536. DOUBLE t; 
  537. {
  538. DOUBLE *p;
  539. int j;
  540.  
  541. hstep = h;
  542. ineq = neq;
  543. /*ccor = hstep * ((DOUBLE *)&precof) [ord];*/
  544. ccor = hstep * ((DOUBLE *)precof) [ord];
  545. jstep = 0;
  546. dsiz = ord + 2;
  547. asiz = neq * dsiz;
  548. order = ord;
  549. ordp1 = ord + 1;
  550.  
  551. setptrs(work);
  552.  
  553. func( t, yn, vn );
  554.  
  555. p = dv;
  556. for( j=0; j<neq; j++ )
  557.     {
  558.     dupdate( p, 0, vn[j] );
  559.     p += dsiz;
  560.     }
  561. }
  562.  
  563.  
  564. /* Calculate and initialize array pointers
  565.  */
  566. setptrs( w )
  567. DOUBLE *w;
  568. {
  569. DOUBLE *p;
  570.  
  571. p = w;
  572. dv = p;
  573. p += asiz;
  574. dvp = p;
  575. p += dsiz;
  576. vp = p;
  577. p += ineq;
  578. yp = p;
  579. p += ineq;
  580. vn = p;
  581. p += ineq;
  582. delta = p;
  583. p += ineq;
  584. sdelta = p;
  585. p += ineq;
  586. adamy0 = p;
  587. }
  588.  
  589.  
  590. /* Subroutines to save and restore the internal
  591.  * state of the integrator.
  592.  */
  593. sinternal(f) /* save state */
  594. register int f; /* file */
  595. {
  596. write( f, &ineq, sizeof(int) );
  597. write( f, &order, sizeof(int) );
  598. write( f, &ordp1, sizeof(int) );
  599. write( f, &asiz, sizeof(int) );
  600. write( f, &dsiz, sizeof(int) );
  601. write( f, &jstep, sizeof(long) );
  602. write( f, &ccor, sizeof(DOUBLE) );
  603. write( f, &hstep, sizeof(DOUBLE) );
  604. }
  605.  
  606. rinternal(f, work) /* restore state */
  607. register int f;
  608. DOUBLE *work;
  609. {
  610. read( f, &ineq, sizeof(int) );
  611. read( f, &order, sizeof(int) );
  612. read( f, &ordp1, sizeof(int) );
  613. read( f, &asiz, sizeof(int) );
  614. read( f, &dsiz, sizeof(int) );
  615. read( f, &jstep, sizeof(long) );
  616. read( f, &ccor, sizeof(DOUBLE) );
  617. read( f, &hstep, sizeof(DOUBLE) );
  618. setptrs( work );
  619. }
  620.  
  621.  
  622. /* Adams-Bashforth-Moulton step.
  623.  */
  624.  
  625. DOUBLE adstep( t, yn, nequat )
  626. DOUBLE *t;
  627. DOUBLE yn[];
  628. int nequat;
  629. {
  630. DOUBLE e, e0, time;
  631. DOUBLE *pdv;
  632. int i, j;
  633. DOUBLE intpol();
  634.  
  635. time = *t;
  636. jstep += 1;
  637.  
  638. /* Do Runge-Kutta for the first ord + 1 steps.
  639.  */
  640. if( jstep <= ordp1 )
  641.     {
  642. /* Take two half steps */
  643.     e0 = Half * hstep;
  644.     rungek(  nequat, time, yn, e0, yn, delta );
  645.     e = time + e0;
  646.     func( e, yn, vn );
  647.     rungek(  nequat, e, yn, e0, yn, delta );
  648.     func( time+hstep, yn, vn );
  649.     pdv = dv;
  650.     for( j=0; j<nequat; j++ )
  651.         {
  652.         dupdate( pdv, (int )jstep, vn[j] );
  653.         pdv += dsiz;
  654.         }
  655.     e = Zero;
  656.     goto done;
  657.     }
  658.  
  659.  
  660.  
  661. /* Predict the next position
  662.  * based on current y' difference table dv[].
  663.  */
  664.     pdv = dv;
  665.     for( i=0; i<nequat; i++ )
  666.         {
  667.         yp[i] = yn[i] + hstep * intpol( pdv, precof, order );
  668.         pdv += dsiz;
  669.         }
  670. /* Evaluate derivatives at the predicted position
  671.  */
  672.     func( time+hstep, yp, vp );
  673.  
  674. /* Correct the predicted position and velocity using the derivatives
  675.  * evalutated at yp.
  676.  */
  677.     pdv = dv;
  678.     for( i=0; i<nequat; i++ )
  679.         {
  680.         vcopy( pdv, dvp, dsiz ); /* y' difference table */
  681.         dupdate( dvp, ordp1, vp[i] );
  682. /* Note, the following line is equivalent to:
  683.  *        yn[i] = yn[i] + hstep * intpol( dap, corcof, ordp1 );
  684.  */
  685.         yn[i] = yp[i] + ccor * dvp[order];
  686.         pdv += dsiz;
  687.         }
  688.  
  689. /* Evaluate derivative at the final position
  690.  */
  691.     func( time+hstep, yn, vn );
  692.  
  693.     e = Zero;
  694.     pdv = dv;
  695.     for( i=0; i<nequat; i++ )
  696.         {
  697.         dupdate( pdv, ordp1, vn[i] );
  698. /* Note, the following line is equivalent to:
  699.  *        yn[i] = yn[i] + hstep * intpol( dap, corcof, ordp1 );
  700.  */
  701.         yn[i] = yp[i] + ccor * pdv[order];
  702. /* Estimate the error
  703.  */
  704. /*        e0 = hstep * pdv[order] * ((DOUBLE *)&corcof) [order];*/
  705.         e0 = hstep * pdv[order] * ((DOUBLE *)corcof) [order];
  706.         if( e0 < 0 )
  707.             e0 = -e0;
  708.         if( e0 > e )
  709.             e = e0;
  710.         pdv += dsiz;
  711.         }
  712. done:
  713. time += hstep;
  714. *t = time;
  715. return( e );
  716. }
  717.